Last updated: 2017-04-02

Readings

  • Garland T Jr, Adolph SC (1994) Why not to do two-species comparative studies: limitations on inferring adaptation. Physiol Zool 67:797-828.

Resources

  • Nunn CL. (2011) The Comparative Approach in Evolutionary Anthropology and Biology. Chicago: University of Chicago Press.
  • Nunn CL, Barton RA (2001) Comparative methods for studying primate adaptation and allometry. Evol Anthropol 10:81–98.

Even more resources

Yet more resources

Non-R resources

R Resources

Trees and tree thinking

  • As soon as you have comparative data, you must consider the relationships among your taxa
    • No longer an option (30+ years of methods)
  • Not concerned in this class with making the tree(s)

You have a tree and data - What do I do now?

Where does your tree come from?

  • Published trees
    • As is
    • Pruned down from published trees
  • Constructed from multiple trees
  • Including tree uncertainty
    • Compare effect of different taxon placements
Trees are hypotheses of relationships. You have hypotheses resting on hypotheses.

Is body size associated with home range area in mammals?

Data from Garland et al. (1992)

  • Body mass (km)
  • Home range area (km2)
  • "49lbr" data

Home range vs. body mass

Home range vs. body mass

Home range vs. body mass

Home range vs. body mass

Why phylogenetic comparative methods?

  1. Meet the assumptions of standard statistical tests
    • Any test you can do using "standard" methods has a phylogenetically informed counterpart.
    • What is your evolutionary model?
  2. Ask and answer nuanced evolutionary questions
    • Trait evolution (continuous or categorical)
    • Trait co-evolution
    • Rates of evolution (identify nodes when rates shift significantly)
    • What is your evolutionary model?

Branch lengths \(\approx\) evolutionary model

"Real" branch lengths:

  • Divergence dates (e.g., fossil calibrated)
  • Genetic or character distance (relative)

Arbitrary branch lengths:

  • All = 1
  • Pagel's ultrametric transformation (1992): total height = deepest nesting

Either can be transformed.

Divergence dated

All = 1

Square-root divergence time

Models of evolution

  • Brownian motion
  • Ornstein-Uhlenbeck (OU)
  • Variable rates (ACDC)
  • Pagel's \(\lambda\)

These models are all accomplished via branch length transformations.

Brownian motion

  • The default, null underlying model for evolution
  • Simplest, most parsimonious model of drift
  • No selection
  • At each step, a trait value is randomly altered based on a normal distribution (mean usually 0) with some variance.

Brownian motion

Simulate 100 Brownian motion random walks with a length of 100 steps.

set.seed(3)
nsim <- 100  # Number of simulations
t <- 0:100   # Time
s2 <- 0.01   # Variance

X <- foreach(i = 1:nsim, .combine = rbind) %do%
  c(0, cumsum(rnorm(n = length(t) - 1, sd = sqrt(s2))))

Brownian motion

Brownian motion

Brownian motion

Brownian motion

Ornstein-Uhlenbeck model

  • Based on Ornstein & Uhlenbeck (1930)
  • Brownian motion doesn't allow for adaptive evolution (drift only)
  • Simplest model for an evolutionary process with selection (deterministic + stochastic)

The OU model describes the motion of a species in the phenotypic space whereby the species moves randomly within the space, but is influenced by a central tendency such that large deviations from the central optimum receive a stronger force back toward the optimum (Blomberg et al. 2003)

Ornstein-Uhlenbeck model (OU)

\[dX(t) = \alpha[\theta - X(t)]dt + \sigma dB(t)\]

  • Change in trait X at time t is a function of selection \(\alpha[\theta - X(t)]dt\) and drift \(\sigma dB(t)\)
    • \(\alpha\) = strength of selection
    • \(\theta\) = optimum
    • \(\sigma\) = magnitude of drift
  • \(\alpha\) and one or more local optima (\(\theta_1\), \(\theta_2\), etc.) can be estimates
  • Models compared with likelihood ratio tests, AIC, etc.

Ornstein-Uhlenbeck model (OU)

Key readings:

  • Martins E, Hansen TF (1996) Phylogenies and the comparative method: A general approach to incorporating phylogenetic information into the analysis of interspecific data. Am Nat 149:646-667.
  • Blomberg SP, Garland T Jr, Ives AR (2003) Testing for phylogenetic signal in comparative data: Behavioral traits are more labile. Evolution 57:717-745.
  • Butler M, King A (2004) Phylogenetic comparative analysis: a modeling approach for adaptive evolution. Am Nat 164:683-695.

OU evolution along a tree

library(phytools)
set.seed(331)
tree <- pbtree(n = 15, scale = 10)

OU evolution along a tree

OU evolution along a tree

\[dX(t) = \alpha[\theta - X(t)]dt + \sigma dB(t)\]

x <- fastBM(tree,
            a = 0,           # Ancestral state at the root node
            theta = 3,       # Optimum
            alpha = 0.2,     # Strength of selection
            s2 = 0.1,      # Variance of the BM process
            internal = TRUE) # Return states for internal nodes

OU evolution along a tree

OU evolution along a tree

OU evolution along a tree

OU evolution along a tree

Estimating an OU model

Rather than simulating data, estimate the values for \(\alpha\) and \(\theta\).

  1. Given a tree and data (body mass in Caribbean anoles)
  2. What model of evolution best fits the data?
    • Brownian motion
    • OU model with a single, clade-wide optimum
    • OU model with three body size optima (S/M/L)
    • OU model with four optima (including an ancestral size)
    • OU model based on multiple colonization events

Estimating an OU model

Evolution toward two local optima:

Estimating an OU model

Estimating an OU model

Estimating an OU model

\(\theta_{small} = 0.25\) mm head length.

ACDC

  • Evolution that either increases or decreases in rate over time.
    • Adaptive radiations followed by phenotypic stasis
    • Proposed by Blomberg et al. (2003)
  • "Early burst" model of Harmon et al. (2010)
  • Didn't really ever catch on as much as OU.

Models of evolution \(\approx\) branch length transformations

Pagel's \(\lambda\): how well does the tree fit the predicted covariance among the trait values

  • \(1\) = Brownian motion
  • \(>1\) = nodes pulled toward the tips (stronger covariance than expected by BM)
  • \(<1\) = nodes pulled toward the root (weaker covariance than expected by BM; more common)

\[\lim_{d\to0} = \mbox{"Star phylogeny"}\]

Pagel's \(\lambda\) branch scaling

Pagel's \(\lambda\) branch scaling

Pagel's \(\lambda\) branch scaling

Pagel's \(\lambda\) branch scaling

Home range vs. body mass

Home range vs. body mass

Working with nlme::gls()

  • We used gls() to fit the squid model, allowing correlated variance on a per-month basis
  • Phylogenetically based correlation variation is a similar issue
  • Choose a phylogenetic correlation structure

Phylogenetic correlation structures

See ?corClasses

  • corBrownian: BM model of correlation
  • corPagel: Pagel's \(\lambda\)
  • corMartins: OU model originally proposed by Martins and Hansen (1997)
  • corBlomberg: ACDC model of Blomberg et al. (2003)

Model 1: Mass, BM

library(nlme)
fm1 <- gls(log_Home_Range ~ log_Mass, data = lbr,
           correlation = corBrownian(phy = tree))
summary(fm1)

Model 1: Mass, BM

## Generalized least squares fit by REML
##   Model: log_Home_Range ~ log_Mass 
##   Data: lbr 
##        AIC      BIC    logLik
##   94.44211 99.99256 -44.22106
## 
## Correlation Structure: corBrownian
##  Formula: ~1 
##  Parameter estimate(s):
## numeric(0)
## 
## Coefficients:
##                 Value Std.Error   t-value p-value
## (Intercept) -1.423845 0.6193642 -2.298882   0.026
## log_Mass     1.261576 0.1767717  7.136753   0.000
## 
##  Correlation: 
##          (Intr)
## log_Mass -0.572
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -1.5273768 -0.3768099  0.2233644  0.5177016  1.5821557 
## 
## Residual standard error: 1.255639 
## Degrees of freedom: 49 total; 47 residual

Model 2: Mass + Clade, BM

fm2 <- gls(log_Home_Range ~ log_Mass + Clade , data = lbr,
           correlation = corBrownian(phy = tree))
summary(fm2)
## Generalized least squares fit by REML
##   Model: log_Home_Range ~ log_Mass + Clade 
##   Data: lbr 
##        AIC      BIC    logLik
##   92.13749 99.45206 -42.06875
## 
## Correlation Structure: corBrownian
##  Formula: ~1 
##  Parameter estimate(s):
## numeric(0)
## 
## Coefficients:
##                      Value Std.Error   t-value p-value
## (Intercept)      -2.233013 0.8012034 -2.787073  0.0077
## log_Mass          1.317072 0.1777494  7.409710  0.0000
## CladeCarnivorans  1.605512 1.0302879  1.558314  0.1260
## 
##  Correlation: 
##                  (Intr) lg_Mss
## log_Mass         -0.557       
## CladeCarnivorans -0.648  0.200
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -1.0443623 -0.1687814  0.1826581  0.6855996  1.0732801 
## 
## Residual standard error: 1.236984 
## Degrees of freedom: 49 total; 46 residual

Model 3: Mass * Clade, BM

fm3 <- gls(log_Home_Range ~ log_Mass * Clade , data = lbr,
           correlation = corBrownian(phy = tree))
summary(fm3)
## Generalized least squares fit by REML
##   Model: log_Home_Range ~ log_Mass * Clade 
##   Data: lbr 
##        AIC      BIC    logLik
##   94.16136 103.1947 -42.08068
## 
## Correlation Structure: corBrownian
##  Formula: ~1 
##  Parameter estimate(s):
## numeric(0)
## 
## Coefficients:
##                               Value Std.Error   t-value p-value
## (Intercept)               -2.322649 0.8699739 -2.669792  0.0105
## log_Mass                   1.352785 0.2200051  6.148882  0.0000
## CladeCarnivorans           1.791226 1.2329695  1.452774  0.1532
## log_Mass:CladeCarnivorans -0.106957 0.3807330 -0.280924  0.7801
## 
##  Correlation: 
##                           (Intr) lg_Mss CldCrn
## log_Mass                  -0.635              
## CladeCarnivorans          -0.706  0.448       
## log_Mass:CladeCarnivorans  0.367 -0.578 -0.536
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -1.0564654 -0.1557668  0.1917349  0.6952910  1.0664477 
## 
## Residual standard error: 1.249557 
## Degrees of freedom: 49 total; 45 residual

Model 4: Mass, OU

Estimating OU parameter alpha often fails (Error in corFactor.corStruct(object))

  • Brute-force optimization of OU alpha parameter
  • Maximize the log-likelihood for a range of models with different values of alpha.
myOU <- function(alpha, dat, tree){
  fm <- gls(log_Home_Range ~ log_Mass, data = dat, 
            correlation = corMartins(alpha, tree, fixed = TRUE))
  return(as.numeric(fm$logLik))
}

alpha <- optimize(f = myOU, interval = c(0, 2), 
                  dat = lbr, tree = tree, maximum = TRUE)
alpha
## $maximum
## [1] 0.01209072
## 
## $objective
## [1] -43.50331

Model 4: Mass, OU

Fit a model with the alpha parameter fixed at the optimized value:

ou_mod <- corMartins(value = alpha$maximum,
                     phy = tree, fixed = TRUE)
fm4 <- gls(log_Home_Range ~ log_Mass, data = lbr,
           correlation = ou_mod)

Model 4: Mass, OU

## Generalized least squares fit by REML
##   Model: log_Home_Range ~ log_Mass 
##   Data: lbr 
##        AIC      BIC    logLik
##   93.00662 98.55706 -43.50331
## 
## Correlation Structure: corMartins
##  Formula: ~1 
##  Parameter estimate(s):
##      alpha 
## 0.01209072 
## 
## Coefficients:
##                 Value Std.Error   t-value p-value
## (Intercept) -1.212025 0.6601609 -1.835954  0.0727
## log_Mass     1.165679 0.1737639  6.708405  0.0000
## 
##  Correlation: 
##          (Intr)
## log_Mass -0.521
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -1.7003954 -0.3764105  0.2171781  0.5543371  1.7814896 
## 
## Residual standard error: 1.066277 
## Degrees of freedom: 49 total; 47 residual

Model comparison with AICc

Models 5 and 6 ANCOVA with OU. Models 7-9, estimate Pagel's \(\lambda\).

Mass Mass + Clade Mass * Clade
Brownian Motion 95.0 93.0 95.6
Stabilizing Selection (OU) 93.5 87.6 90.3
Pagel's \(\lambda\) 96.7 93.5 96.2

Best performing model

OU model:

\[Home~Range \sim Mass + Clade\]

## Generalized least squares fit by REML
##   Model: log_Home_Range ~ log_Mass + Clade 
##   Data: lbr 
##        AIC      BIC    logLik
##   86.67971 93.99427 -39.33985
## 
## Correlation Structure: corMartins
##  Formula: ~1 
##  Parameter estimate(s):
##      alpha 
## 0.03533854 
## 
## Coefficients:
##                      Value Std.Error   t-value p-value
## (Intercept)      -1.919568 0.4851383 -3.956743  0.0003
## log_Mass          1.226691 0.1686608  7.273123  0.0000
## CladeCarnivorans  1.391452 0.4100937  3.393010  0.0014
## 
##  Correlation: 
##                  (Intr) lg_Mss
## log_Mass         -0.855       
## CladeCarnivorans -0.707  0.475
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -1.8909506 -0.3117712  0.2184212  1.0933805  1.7774159 
## 
## Residual standard error: 0.6911628 
## Degrees of freedom: 49 total; 46 residual

No quiz

Lecture 11-3

References

Blomberg, S. P., T. Garland Jr, and A. R. Ives. 2003. Testing for Phylogenetic Signal in Comparative Data: Behavioral Traits Are More Labile. Evolution 57:717–745.

Butler, M., and A. King. 2004. Phylogenetic Comparative Analysis: A Modeling Approach for Adaptive Evolution. American Naturalist 164:683–695.

Garland, T., Jr, and S. C. Adolph. 1994. Why Not to Do Two-Species Comparative Studies: Limitations on Inferring Adaptation. Physiological Zoology 67:797–828.

Garland, T., Jr, P. H. Harvey, and A. R. Ives. 1992. Procedures for the Analysis of Comparative Data Using Phylogenetically Independent Contrasts. Systematic Biology 41:18–32.

Hansen, T. F., and E. Martins. 1996. Translating Between Microevolutionary Process and Macroevolutionary Patterns: The Correlation Structure of Interspecific Data. Evolution 50:1404–1417.

Harmon, L. J., J. B. Losos, T. J. Davies, R. G. Gillespie, J. L. Gittleman, W. Bryan Jennings, K. H. Kozak, M. A. McPeek, F. Moreno-Roark, T. J. Near, A. Purvis, R. E. Ricklefs, D. Schluter, J. A. Schulte Ii, O. Seehausen, B. L. Sidlauskas, O. Torres-Carvajal, J. T. Weir, and A. Ø. Mooers. 2010. Early Bursts of Body Size and Shape Evolution Are Rare in Comparative Data. Evolution 64:2385–2396.

Martins, E., and T. F. Hansen. 1997. Phylogenies and the Comparative Method: A General Approach to Incorporating Phylogenetic Information into the Analysis of Interspecific Data. American Naturalist 149:646–667.

Nunn, C. L. 2011. The Comparative Approach in Evolutionary Anthropology and Biology. University of Chicago Press, Chicago.

Nunn, C. L., and R. A. Barton. 2001. Comparative Methods for Studying Primate Adaptation and Allometry. Evolutionary Anthropology 10:81–98.

Uhlenbeck, G. E., and L. S. Ornstein. 1930. On the Theory of the Brownian Motion. Physical Review 36:823–841.